#include "numuCCQESelection.hxx"
#include "numuCCSelection.hxx"
#include "baseSelection.hxx"
#include "CutUtils.hxx"

//reduced FGD1 FV 
const   Float_t _FVdefminFGD1[3] = {57.66, 57.66, 20.925};
const   Float_t _FVdefmaxFGD1[3] = {57.66, 57.66, 20.925};

//BarrelECal FV for stopping point cut 
const   Float_t _FVdefminTECALLeft[3] = {0.0, 40.0, 40.0};
const   Float_t _FVdefmaxTECALLeft[3] = {10.0, 40.0, 40.0};

const   Float_t _FVdefminTECALRight[3] = {10, 40.0, 40.0};
const   Float_t _FVdefmaxTECALRight[3] = {0.0, 40.0, 40.0};

const   Float_t _FVdefminTECALTop[3] = {40.0, 0.0, 40.0};
const   Float_t _FVdefmaxTECALTop[3] = {40.0, 10.0, 40.0};

const   Float_t _FVdefminTECALBottom[3] = {40.0, 10.0, 40.0};
const   Float_t _FVdefmaxTECALBottom[3] = {40.0, 0.0, 40.0};

//SMRD FV for stopping point cut 
const   Float_t _FVdefminSMRDLeft[3] = {0.0, 0.0, 125.0};
const   Float_t _FVdefmaxSMRDLeft[3] = {60.0, 0.0, 125.0};

const   Float_t _FVdefminSMRDRight[3] = {60.0, 0.0, 125.0};
const   Float_t _FVdefmaxSMRDRight[3] = {0.0, 0.0, 125.0};

const   Float_t _FVdefminSMRDTop[3] = {0.0, 0.0, 125.0};
const   Float_t _FVdefmaxSMRDTop[3] = {0.0, 60.0, 125.0};

const   Float_t _FVdefminSMRDBottom[3] = {0.0, 60.0, 125.0};
const   Float_t _FVdefmaxSMRDBottom[3] = {0.0, 0.0, 125.0};


//********************************************************************
numuCCQESelection::numuCCQESelection(bool forceBreak): SelectionBase(forceBreak) {
//********************************************************************

    // define the steps
    DefineSteps();  
}

//********************************************************************
void numuCCQESelection::DefineSteps(){
//********************************************************************

    // Cuts must be added in the right order
    // last "true" means the step sequence is broken if cut is not passed (default is "false")
    // Add inclusive CC steps to all branches

    AddStep(StepBase::kCut,    "event quality",        new EventQualityCut(),           true);
    AddStep(StepBase::kCut,    "> 0 tracks ",          new TotalMultiplicityCut(),      true);  
    AddStep(StepBase::kAction, "find leading tracks",  new FindLeadingTracksAction());  
    AddStep(StepBase::kAction, "find tracker tracks",  new FindTrackerTracksAction());
    AddStep(StepBase::kAction, "find candidates",      new FindCandidatesAction());
    AddStep(StepBase::kAction, "fill summary",         new FillSummaryAction_numuCCQE());

    // Add a split to the trunk with 2 branches 
    AddSplit(2);

    // First branch is for cases 0,1,2
    AddStep(0,StepBase::kAction, "find vertex",        new FindVertexAction());  
    AddStep(0,StepBase::kCut,    "quality+fiducial",   new TrackQualityFiducialCut(),   true);  
    AddStep(0,StepBase::kAction, "find veto track",    new FindVetoTrackAction());
    AddStep(0,StepBase::kCut,    "veto",               new ExternalVetoCut());
    AddStep(0,StepBase::kAction, "find oofv track",    new FindOOFVTrackAction());
    // AddStep(0,StepBase::kCut,    "External FGD1",      new ExternalFGD1lastlayersCut());
    AddStep(0,StepBase::kCut,    "muon PID",           new MuonPIDCut());
    AddStep(0,StepBase::kCut,    "1 neg track",        new OneNegativeTrackCut());
    AddStep(0,StepBase::kAction, "find michel e-",     new FindMichelElectronsAction());
    AddStep(0,StepBase::kCut,    "no michel ele",      new MichelElectronCut());
    AddStep(0,StepBase::kAction, "find tracker tracks",new FindTrackerTracksAction());


    // Add a split to branch 0 with 2 sub-branches 
    AddSplit(2,0);

    //--------- Case 0 ------------ 
    AddStep(0,0, StepBase::kCut, "one tracker track",    new OneTrackerTrackCut());
    //  DisableStep(0,"ExternalFGD1");



    // This is for cases 1 and 2

    AddStep(0,1,StepBase::kCut,    "two tracker tracks",       new TwoTrackerTracksCut());
    //if (_findCorrectSense)
    AddStep(0,1,StepBase::kAction, "find correct sense",       new FindCorrectProtonSenseAction());


    // Add a split to branch (0,1) with 2 sub-branches 
    AddSplit(2,0,1);

    //--------- Case 1 ------------ 

    AddStep(0,1,0,StepBase::kCut, "1 pos track",        new OnePositiveTrackCut());
    AddStep(0,1,0,StepBase::kCut, "common vertex",      new CommonVertexCut1());
    AddStep(0,1,0,StepBase::kCut, "proton PID",         new ProtonPIDCut());
    //  DisableStep(1,0,"ExternalFGD1");

    //--------- Case 2 ------------ 
    AddStep(0,1,1,StepBase::kCut,       "1 contained FGD iso", new OneFGDContainedTrackCut());
    AddStep(0,1,1,StepBase::kCut,       "common vertex",       new CommonVertexCut2());
    AddStep(0,1,1,StepBase::kCut,       "FGD proton PID",      new FGDProtonPIDCut());
    AddStep(0,1,1,StepBase::kCut,       "stopping proton",     new StoppingProtonCut());
    //  DisableStep(2,"ExternalFGD1");


    //--------- Case 3 ------------ 

    AddStep(1, StepBase::kAction, "find proton vertex",    new FindProtonVertexAction());  // Find the vertex (HMP track for the moment)
    AddStep(1, StepBase::kCut,    "quality+fiducial",      new TrackQualityFiducialCut3(), true);
    AddStep(1, StepBase::kAction, "find veto track",       new FindVetoTrackAction3());
    AddStep(1, StepBase::kCut,    "veto",                  new ExternalVetoCut3());
    AddStep(1, StepBase::kCut,    "proton PID",            new ProtonPIDCut());
    AddStep(1, StepBase::kCut,    "1 pos track",           new OnePositiveTrackCut());
    AddStep(1, StepBase::kCut,    "two tracker tracks",    new TwoTrackerTracksCut());

    //  if (_findCorrectSense)
    AddStep(1, StepBase::kAction, "find correct muon sense", new FindCorrectMuonSenseAction());
    AddStep(1, StepBase::kCut,    "1 FGD1 iso",              new OneFGD1TrackCut()); // in FGD1 volume (change it!!!)
    AddStep(1, StepBase::kCut,    "long FGD track",          new LongFGDTrackCut());
    AddStep(1, StepBase::kCut,    "FGD muon PID",            new FGDMuonPIDCut());
    AddStep(1, StepBase::kCut,    "common vertex",           new CommonVertexCut3());
    AddStep(1, StepBase::kCut,    "stopping muon",           new StoppingMuonCut());



    SetBranchAlias(0,"muTPC",0,0);
    SetBranchAlias(1,"muTPC-pTPC",0,1,0);
    SetBranchAlias(2,"muTPC-pFGD",0,1,1);
    SetBranchAlias(3,"muFGD-pTPC",1);
}

//********************************************************************
bool numuCCQESelection::FillEventSummary(const AnaBoxB& box, AnaEventB& event, Int_t allCutsPassed[]){
//********************************************************************

    (void)box;

    // Should set event sample enumeration here
    if(allCutsPassed[0] || allCutsPassed[1] || allCutsPassed[2] || allCutsPassed[3]) event.Summary->EventSample = nd280Samples::kFGD1NuMuCCQE;
    return (event.Summary->EventSample != nd280Samples::kUnassigned);
}


//********************************************************************
bool FillSummaryAction_numuCCQE::Apply(const AnaEventB& event, AnaBoxB& box) const{
//********************************************************************

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if(!ccqebox->MuonCandidate) return 1;

    // Have the muon candidate, so fill lepton candidate summary info
    event.Summary->LeptonCandidate[nd280Samples::kFGD1NuMuCCQE] = ccqebox->MuonCandidate;

    // We have the HM negative and positive candidates now, so can fill the summary information
    if(box.HMNtrack){
        for(int i = 0; i < 4; ++i){
            event.Summary->VertexPosition[nd280Samples::kFGD1NuMuCCQE][i] = box.HMNtrack->PositionStart[i];
        }
        if(box.HMNtrack->TrueTrack) event.Summary->TrueVertex[nd280Samples::kFGD1NuMuCCQE] = box.HMNtrack->TrueTrack->TrueVertex;
    }
    else if(box.HMPtrack){
        for(int i = 0; i < 4; ++i){
            event.Summary->VertexPosition[nd280Samples::kFGD1NuMuCCQE][i] = box.HMPtrack->PositionStart[i];
        }

        if(box.HMPtrack->TrueTrack) event.Summary->TrueVertex[nd280Samples::kFGD1NuMuCCQE] = box.HMPtrack->TrueTrack->TrueVertex;
    }
    return 1;
}

//**************************************************
bool OneNegativeTrackCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    // Exactly one negative track with good quality and in FV
    if (box.nNegativeTPCtracks==1) 
        return true;
    else 
        return false;
}

//**************************************************
bool MichelElectronCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    // Temporarily, a Float_t cut to make sure it works with OldFlatTree, FlatTree and oaAnalysis
    if (box.nFGD1MichelElectrons == 0)
        return true;
    else
        return false;
}

//**************************************************
bool OneTrackerTrackCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (ccqebox->TrackerTracks.size() == 1) 
        return true;
    else 
        return false;
}

//**************************************************
bool TwoTrackerTracksCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (ccqebox->TrackerTracks.size() == 2) 
        return true;
    else 
        return false;
}

//**************************************************
bool OnePositiveTrackCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    // Exactly one positive track with good quality and in FV
    if (box.nPositiveTPCtracks == 1) 
        return true;
    else 
        return false;
}

//**************************************************
bool CommonVertexCut1::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (!ccqebox->MuonCandidate || !ccqebox->ProtonCandidate) return false;

    return cutUtils::CommonVertexCut(*(ccqebox->MuonCandidate), *(ccqebox->ProtonCandidate), 
            _cutX, _cutY, _cutZ);
}

//**************************************************
bool CommonVertexCut2::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (!ccqebox->MuonCandidate || !ccqebox->ProtonCandidate) return false;

    return cutUtils::CommonVertexCut(*(ccqebox->MuonCandidate), *(ccqebox->ProtonCandidate), 
            _cutX, _cutY, _cutZ);
}

//**************************************************
bool CommonVertexCut3::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (!ccqebox->MuonCandidate || !ccqebox->ProtonCandidate) return false;

    return cutUtils::CommonVertexCut(*(ccqebox->MuonCandidate), *(ccqebox->ProtonCandidate), 
            _cutX, _cutY, _cutZ);
}

//**************************************************
bool ProtonPIDCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    if (!box.HMPtrack) return false;

    // use different cut value w.r.t. to cutUtils::ProtonPIDCut()  
    if (anaUtils::GetPIDLikelihood(*(box.HMPtrack),2)>0.5)
        return true;
    else
        return false;

}

//**************************************************
bool OneFGDContainedTrackCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    // Only one FGD-iso track (with no other detectors) in FGD1 Volume
    if (ccqebox->FGDtracksCrossingFGD1.size() != 1) return false;

    if (anaUtils::InDetVolume(SubDetId::kFGD1,ccqebox->FGDtracksCrossingFGD1[0]->PositionStart) && 
            anaUtils::TrackUsesOnlyDet(*ccqebox->FGDtracksCrossingFGD1[0], SubDetId::kFGD1))
        return true;
    else
        return false;
}

//**************************************************
bool FGDProtonPIDCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (!ccqebox->ProtonCandidate) return false;
    return std_numuccqe_actions::FGDProtonPIDCut(*(ccqebox->ProtonCandidate));
}


//********************************************************************
bool ExternalVetoCut3::Apply(const AnaEventB& event, AnaBoxB& box) const{
//********************************************************************

    (void)event;

    AnaBoxCCQE* ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (!box.HMPtrack) return false;
    return cutUtils::ExternalVetoCut(*(box.HMPtrack), ccqebox->VetoTrack3);
}

//**************************************************
bool OneFGD1TrackCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    // Only one iso track crossing FGD1. We are implicitly requiring no FGD2 iso-tracks since we ask for two tracker tracks
    if (ccqebox->FGDtracksCrossingFGD1.size() != 1) return false;

    if (anaUtils::InDetVolume(SubDetId::kFGD1,ccqebox->FGDtracksCrossingFGD1[0]->PositionStart))
        return true;
    else
        return false;

}

//**************************************************
bool LongFGDTrackCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (!ccqebox->MuonCandidate) return false;

    Float_t length = 0;

    if(ccqebox->MuonCandidate->nFGDSegments>0) 
        length = ccqebox->MuonCandidate->FGDSegments[0]->X;

    Float_t cut=500;

    if (!anaUtils::TrackUsesOnlyDet(*(ccqebox->MuonCandidate), SubDetId::kFGD1) || length>cut)
        return true;
    else 
        return false;
}

//**************************************************
bool FGDMuonPIDCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (!ccqebox->MuonCandidate) return false;

    // apply it only if the FGD track is contained 
    if (!anaUtils::TrackUsesOnlyDet(*(ccqebox->MuonCandidate), SubDetId::kFGD1)) return true;

    return std_numuccqe_actions::FGDMuonPIDCut(*(ccqebox->MuonCandidate));
}

//**************************************************
bool StoppingProtonCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    // Needed to compute the momentum by range

    if(!ccqebox->ProtonCandidate) return false;

    // apply it only if the FGD track is contained 
    if (!anaUtils::TrackUsesOnlyDet(*(ccqebox->ProtonCandidate), SubDetId::kFGD1)) return false;

    if (std_numuccqe_actions::TrackEndsInFgd1FV(*(ccqebox->ProtonCandidate)))
        return true;
    else 
        return false;
}

//**************************************************
bool TrackQualityFiducialCut3::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************
    // Needed to compute the momentum by range

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (ccqebox->ProtonVertex && cutUtils::FiducialCut(*(box.HMPtrack), SubDetId::kFGD1)) 
        //if (ccqebox->ProtonVertex) // to remove the fiducial volume cut 
        return true;
    else 
        return false;
}



//**************************************************
bool StoppingMuonCut::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************
    // Needed to compute the momentum by range

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (!ccqebox->MuonCandidate) return false;

    // if the FGD track is only in FGD1 it should be contained
    if (anaUtils::TrackUsesOnlyDet(*(ccqebox->MuonCandidate), SubDetId::kFGD1)){
        if (std_numuccqe_actions::TrackEndsInFgd1FV(*(ccqebox->MuonCandidate)))
            return true;
        else 
            return false;
    }
    else{
        // if the FGD track has other detectors it should be contained in those
        if (std_numuccqe_actions::TrackEndsInBarrelECalFV(*(ccqebox->MuonCandidate)))
            return true;
        else if (std_numuccqe_actions::TrackEndsInSMRDFV(*(ccqebox->MuonCandidate)))
            return true;
        else
            return false;
    }
}

//**************************************************
bool FindMichelElectronsAction::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    cutUtils::FindFGDMichelElectrons(event,box);

    return true;
}

//**************************************************
bool FindCandidatesAction::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (box.HMNtrack)
        ccqebox->MuonCandidate = box.HMNtrack;
    else if (ccqebox->FGDtracksCrossingFGD1.size()>0)
        ccqebox->MuonCandidate = ccqebox->FGDtracksCrossingFGD1[0]; //tracks crossing FGD1

    if (box.HMPtrack)
        ccqebox->ProtonCandidate = box.HMPtrack;
    else if (box.HMNtrack && ccqebox->FGDtracksCrossingFGD1.size()>0)
        ccqebox->ProtonCandidate = ccqebox->FGDtracksCrossingFGD1[0];  //tracks crossing FGD1

    return true;
}

//**************************************************
bool FindProtonVertexAction::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    // reset the vertex 
    ccqebox->ProtonVertex = NULL;

    if (!box.HMPtrack) return false;

    ccqebox->ProtonVertex = new AnaVertexB();
    utils::CreateArray(ccqebox->ProtonVertex->Tracks, 1);
    ccqebox->ProtonVertex->Tracks[ccqebox->ProtonVertex->nTracks] = box.HMPtrack;
    ccqebox->ProtonVertex->nTracks++;

    cutUtils::AssignArray(box.HMPtrack->PositionStart, ccqebox->ProtonVertex->Position, 4);

    if (box.HMPtrack->TrueTrack) {
        ccqebox->ProtonVertex->TrueVertex = box.HMPtrack->TrueTrack->TrueVertex;
    }

    return true;
}

//**************************************************
bool FindVetoTrackAction3::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    ccqebox->VetoTrack3 = cutUtils::FindVetoTrack(event, *(box.HMPtrack));  
    return true;
}

//**************************************************
bool FindTrackerTracksAction::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    // all Tracker tracks
    ccqebox->TrackerTracks.clear();

    std_numuccqe_actions::FindTracksWithTPC(event,*ccqebox, true, SubDetId::kFGD1);
    std_numuccqe_actions::FindTracksWithFGDAndNoTPC(event, *ccqebox, SubDetId::kFGD1);

    return true;
}

//**************************************************
bool FindCorrectProtonSenseAction::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (ccqebox->ProtonCandidate)
        // it must reverse the track only if it is a FGD track (no TPC segments)
        if (!anaUtils::TrackUsesDet(*(ccqebox->ProtonCandidate), SubDetId::kTPC))
            std_numuccqe_actions::FindCorrectSense(ccqebox->ProtonCandidate, box.Vertex);

    return true;
}

//**************************************************
bool RevertToOriginalProtonSenseAction::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (ccqebox->ProtonCandidate)
        std_numuccqe_actions::RevertToOriginalSense(ccqebox->ProtonCandidate);

    return true;
}

//**************************************************
bool FindCorrectMuonSenseAction::Apply(const AnaEventB& event, AnaBoxB& box) const{
//**************************************************

    (void)event;

    AnaBoxCCQE * ccqebox = static_cast<AnaBoxCCQE*>(&box); 

    if (ccqebox->MuonCandidate)
        // it must reverse the track only if it is a FGD track (no TPC segments)
        if (!anaUtils::TrackUsesDet(*(ccqebox->MuonCandidate), SubDetId::kTPC))
            std_numuccqe_actions::FindCorrectSense(ccqebox->MuonCandidate, ccqebox->ProtonVertex);

    return true;
}

//**************************************************
bool std_numuccqe_actions::FGDProtonPIDCut(const AnaTrackB& track){
//**************************************************

    if (track.nFGDSegments == 0) return false;

    Float_t cut = -4;

    if (track.FGDSegments[0]->Pullp > cut)
        return true;
    else    
        return false;
}

//**************************************************
bool std_numuccqe_actions::FGDMuonPIDCut(const AnaTrackB& track){
//**************************************************

    if (track.nFGDSegments == 0) return false;

    Float_t cut = -5;

    if (track.FGDSegments[0]->Pullmu > cut)
        return true;
    else    
        return false;
}

//**************************************************
bool std_numuccqe_actions::FGDTrackQualityCut(const AnaTrackB& track){
//**************************************************

    if (track.nFGDSegments == 0) return false;

    //  Float_t cut = 3;

    //  if (track.FGDSegments[0]->NHits > cut)
    return true;
    //  else    
    //    return false;
}

//********************************************************************
void std_numuccqe_actions::FindCorrectSense(AnaTrackB* track, AnaVertexB* vertex) {
//********************************************************************

    if (!track) return;
    if (!vertex) return;

    Float_t distStart = cutUtils::GetSeparationSquared(track->PositionStart, vertex->Position);
    Float_t distEnd   = cutUtils::GetSeparationSquared(track->PositionEnd, vertex->Position);

    if (distEnd < distStart){
        // reverse the track
        Float_t posTemp[4] = {track->PositionStart[0], track->PositionStart[1], track->PositionStart[2], track->PositionStart[3]};
        //    TVector3       dirTemp = track->DirectionStart;
        cutUtils::AssignArray(track->PositionEnd, track->PositionStart, 4);
        cutUtils::AssignArray(posTemp, track->PositionEnd, 4);
        //    track->DirectionStart   = - track->DirectionEnd;
        //    track->DirectionEnd     = - dirTemp;    
        track->Charge           *= -1;    
    }

}

//********************************************************************
void std_numuccqe_actions::RevertToOriginalSense(AnaTrackB* track) {
//********************************************************************

    cutUtils::AssignArray(track->Original->PositionStart, track->PositionStart, 4);
    cutUtils::AssignArray(track->Original->PositionEnd, track->PositionEnd, 4);
    cutUtils::AssignArray(track->Original->DirectionStart, track->DirectionStart, 3);
    track->Charge           = track->Original->Charge;
}

//********************************************************************
void std_numuccqe_actions::FindTracksWithTPC(const AnaEventB& event, AnaBoxCCQE& ccqebox, bool qualityCut, SubDetId::SubDetEnum det){
//********************************************************************

    std::vector<AnaTrackB*> allTracks[5];

    // all tracks using TPC
    AnaTrackB* selTracks[100];
    int nTPC = anaUtils::GetAllTracksUsingTPC(event, selTracks);

    //loop over tpc tracks
    for (Int_t i=0;i<nTPC; ++i){
        AnaTrackB* track = selTracks[i];

        // all TPC tracks
        allTracks[4].push_back(track);

        // all Tracker tracks
        ccqebox.TrackerTracks.push_back(track);

        //check if the track starts in FGD1 volume
        if (anaUtils::InDetVolume(SubDetId::kFGD1,track->PositionStart)) {
            if (qualityCut){
                //apply the quality cut: good quality tracks
                if (cutUtils::TrackQualityCut(*track)) { 
                    // Apply the fiducial cut: save tracks out FV
                    if (!cutUtils::FiducialCut(*track, det)) {

                        if (track->Charge == -1)
                            allTracks[0].push_back(track);
                        else
                            allTracks[1].push_back(track);
                    }
                }
                //bad quality tracks
                else {
                    // Apply the fiducial cut: save tracks in FV
                    if (cutUtils::FiducialCut(*track, det))
                        allTracks[2].push_back(track);
                    // save tracks out FV
                    else
                        allTracks[3].push_back(track);
                }
            }
        }
    }

    // Fill the box
    ccqebox.NegativeTPCtracksOutFV = allTracks[0];
    ccqebox.PositiveTPCtracksOutFV = allTracks[1];
    ccqebox.BadQuaTPCtracksInFV = allTracks[2];
    ccqebox.BadQuaTPCtracksOutFV = allTracks[3];
    ccqebox.TPCtracks = allTracks[4];
}

//********************************************************************
void std_numuccqe_actions::FindTracksWithFGDAndNoTPC(const AnaEventB& event, AnaBoxCCQE& ccqebox, SubDetId::SubDetEnum det){
//********************************************************************

    std::vector<AnaTrackB*> allTracks[7];

    // Get all tracks with FGD and no TPC
    AnaTrackB* allFgdTracks[100];
    int nFgdTracks= anaUtils::GetAllTracksUsingFGDAndNoTPC(event, allFgdTracks);

    //loop over tracks
    for (int i=0;i<nFgdTracks;i++){
        AnaTrackB* track = allFgdTracks[i];

        // All tracks using FGD and no TPC
        allTracks[6].push_back(track);      

        // All Tracker tracks
        ccqebox.TrackerTracks.push_back(track);

        // Tracks using FGD1 (and no TPC)
        if (anaUtils::TrackUsesDet(*track,SubDetId::kFGD1)) 
            allTracks[5].push_back(track);      

        //check if the track starts in FGD1 volume
        if (anaUtils::InDetVolume(SubDetId::kFGD1,track->PositionStart)) {
            allTracks[0].push_back(track);      

            // in FV tracks
            if (cutUtils::FiducialCut(*track, det)) {
                // contained FGD tracks
                if (anaUtils::TrackUsesOnlyDet(*track, SubDetId::kFGD1))
                    allTracks[1].push_back(track);
                // not-contained FGD tracks
                else
                    allTracks[2].push_back(track);      
            }
            // out FV tracks
            else {
                // contained FGD tracks
                if (anaUtils::TrackUsesOnlyDet(*track, SubDetId::kFGD1))
                    allTracks[3].push_back(track);
                // not-contained FGD tracks
                else 
                    allTracks[4].push_back(track);      
            }
        }
    }

    // Fill the numuccqe_box
    ccqebox.FGDtracksInFGD1Volume = allTracks[0];
    ccqebox.FGDContTracksInFV = allTracks[1];
    ccqebox.FGDNotContTracksInFV = allTracks[2];
    ccqebox.FGDContTracksOutFV = allTracks[3];
    ccqebox.FGDNotContTracksOutFV = allTracks[4];
    ccqebox.FGDtracksCrossingFGD1 = allTracks[5];
    ccqebox.FGDtracks = allTracks[6];
}

//**************************************************
bool std_numuccqe_actions::TrackEndsInBarrelECal(const AnaTrackB& track) {
//**************************************************

    if(anaUtils::InDetVolume(SubDetId::kTopTECAL,    track.PositionEnd))
        return true;

    if(anaUtils::InDetVolume(SubDetId::kBottomTECAL, track.PositionEnd))
        return true;

    if(anaUtils::InDetVolume(SubDetId::kLeftTECAL,   track.PositionEnd))
        return true;

    if(anaUtils::InDetVolume(SubDetId::kRightTECAL,  track.PositionEnd))
        return true;

    return false;
}

//**************************************************
bool std_numuccqe_actions::TrackEndsInSMRD(const AnaTrackB& track) {
//**************************************************

    if (anaUtils::InDetVolume(SubDetId::kTopSMRD,    track.PositionEnd))
        return true;

    if (anaUtils::InDetVolume(SubDetId::kBottomSMRD, track.PositionEnd))
        return true;

    if (anaUtils::InDetVolume(SubDetId::kLeftSMRD,   track.PositionEnd))
        return true;

    if (anaUtils::InDetVolume(SubDetId::kRightSMRD,  track.PositionEnd))
        return true;

    return false;
}


//**************************************************
bool std_numuccqe_actions::TrackEndsInBarrelECalFV(const AnaTrackB& track){
//**************************************************

    if(anaUtils::InFiducialVolume(SubDetId::kTopTECAL, track.PositionEnd,    _FVdefminTECALTop,    _FVdefmaxTECALTop))
        return true;

    if(anaUtils::InFiducialVolume(SubDetId::kBottomTECAL, track.PositionEnd, _FVdefminTECALBottom, _FVdefmaxTECALBottom))
        return true;

    if(anaUtils::InFiducialVolume(SubDetId::kLeftTECAL, track.PositionEnd,   _FVdefminTECALLeft,   _FVdefmaxTECALLeft))
        return true;

    if(anaUtils::InFiducialVolume(SubDetId::kRightTECAL, track.PositionEnd,  _FVdefminTECALRight,  _FVdefmaxTECALRight))
        return true;

    return false;
}

//**************************************************
bool std_numuccqe_actions::TrackEndsInSMRDFV(const AnaTrackB& track) {
//**************************************************

    if(anaUtils::InFiducialVolume(SubDetId::kTopSMRD, track.PositionEnd,    _FVdefminSMRDTop,    _FVdefmaxSMRDTop))
        return true;

    if(anaUtils::InFiducialVolume(SubDetId::kBottomSMRD, track.PositionEnd, _FVdefminSMRDBottom, _FVdefmaxSMRDBottom))
        return true;

    if(anaUtils::InFiducialVolume(SubDetId::kLeftSMRD, track.PositionEnd,   _FVdefminSMRDLeft,   _FVdefmaxSMRDLeft))
        return true;

    if(anaUtils::InFiducialVolume(SubDetId::kRightSMRD, track.PositionEnd,  _FVdefminSMRDRight,  _FVdefmaxSMRDRight))
        return true;

    return false;
}

//**************************************************
bool std_numuccqe_actions::TrackEndsInFgd1FV(const AnaTrackB& track) {
//**************************************************

    return anaUtils::InFiducialVolume(SubDetId::kFGD1, track.PositionEnd, _FVdefminFGD1, _FVdefmaxFGD1);
}
